diff --git a/QP_solver/include/CGAL/QP_solver.h b/QP_solver/include/CGAL/QP_solver.h
index 273346e4c0e..e77661f10cf 100644
--- a/QP_solver/include/CGAL/QP_solver.h
+++ b/QP_solver/include/CGAL/QP_solver.h
@@ -1340,7 +1340,8 @@ private:
   // validity checks:
   bool is_solution_feasible_for_auxiliary_problem() const;
   bool is_solution_optimal_for_auxiliary_problem() const;
-  bool is_solution_feasible() const;
+  bool is_solution_feasible() const;  
+  bool is_solution_infeasible() const;
   bool is_solution_optimal() const;
   bool is_value_correct() const;
   bool is_solution_unbounded() const;
diff --git a/QP_solver/include/CGAL/QP_solver/Initialization.h b/QP_solver/include/CGAL/QP_solver/Initialization.h
index dd43a3aeadd..fb5ac2cac14 100644
--- a/QP_solver/include/CGAL/QP_solver/Initialization.h
+++ b/QP_solver/include/CGAL/QP_solver/Initialization.h
@@ -195,7 +195,7 @@ init_x_O_v_i()
   // and so we initialize them to zero (if the bound on the variable
   // allows it), or to the variable's lower or upper bound:
   for (int i = 0; i < qp_n; ++i) {
-    CGAL_qpe_assertion(qp_l[i]<=qp_u[i] || !*(qp_fl+i) || !*(qp_fu+i));
+    CGAL_qpe_assertion( !*(qp_fl+i) || !*(qp_fu+i) || qp_l[i]<=qp_u[i]);
 
     if (*(qp_fl+i))                    // finite lower bound?
       if (*(qp_fu+i))                  // finite lower and finite upper bound?
diff --git a/QP_solver/include/CGAL/QP_solver/QP_functions_impl.h b/QP_solver/include/CGAL/QP_solver/QP_functions_impl.h
index e74fa8e9dbe..0bbad36bee5 100644
--- a/QP_solver/include/CGAL/QP_solver/QP_functions_impl.h
+++ b/QP_solver/include/CGAL/QP_solver/QP_functions_impl.h
@@ -135,7 +135,8 @@ namespace QP_functions_detail {
     int n = p.n();
     bool empty_D = true;
     for (int i=0; i<n; ++i, ++it) {
-      for (int j=0; j<n; ++j)
+      // handle only entries on/below diagonal
+      for (int j=0; j<i+1; ++j)
 	if (!CGAL::is_zero(*(*it + j))) {
 	  if (empty_D) {
 	    // first time we see a nonzero entry
@@ -143,6 +144,9 @@ namespace QP_functions_detail {
 	    empty_D = false;
 	  }
 	  out << "  x" << i << "  x" << j << "  " << *(*it + j) << "\n";
+	  // QMATRIX format prescribes symmetric matrix
+	  if (i != j)
+	    out << "  x" << j << "  x" << i << "  " << *(*it + j) << "\n";
 	}
     }
   }
@@ -198,42 +202,41 @@ namespace QP_functions_detail {
 		  << *r1 << " vs. " <<  *r2 << std::endl;	
 	return_val = false;
       }
-    // check fl
+    // check fl, l
     typename Quadratic_program1::FL_iterator fl1 = qp1.fl();
     typename Quadratic_program2::FL_iterator fl2 = qp2.fl();
-    for (int j=0; j<n; ++j, ++fl1, ++fl2)
+    typename Quadratic_program1::L_iterator l1 = qp1.l();
+    typename Quadratic_program2::L_iterator l2 = qp2.l();
+    for (int j=0; j<n; ++j, ++fl1, ++fl2, ++l1, ++l2) {
       if (*fl1 != *fl2) {
 	std::cerr << "Equality test fails with fl[" << j << "]: "
 		  << *fl1 << " vs. " <<  *fl2 << std::endl;	
 	return_val = false;
       }
-    // check l
-    typename Quadratic_program1::L_iterator l1 = qp1.l();
-    typename Quadratic_program2::L_iterator l2 = qp2.l();
-    for (int j=0; j<n; ++j, ++l1, ++l2)
-      if (*l1 != *l2) {
+      if ((*fl1 == true) && (*l1 != *l2)) {
 	std::cerr << "Equality test fails with l[" << j << "]: "
 		  << *l1 << " vs. " <<  *l2 << std::endl;
 	return_val = false;
       }
-    // check fu
+    }
+    
+    // check fu, u 
     typename Quadratic_program1::FU_iterator fu1 = qp1.fu();
     typename Quadratic_program2::FU_iterator fu2 = qp2.fu();
-    for (int j=0; j<n; ++j, ++fu1, ++fu2)
+    typename Quadratic_program1::U_iterator u1 = qp1.u();
+    typename Quadratic_program2::U_iterator u2 = qp2.u();
+    for (int j=0; j<n; ++j, ++fu1, ++fu2, ++u1, ++u2) {
       if (*fu1 != *fu2) {
 	std::cerr << "Equality test fails with fu[" << j << "]: "
 		  << *fu1 << " vs. " <<  *fu2 << std::endl;
 	return_val = false;
       }
-    // check u
-    typename Quadratic_program1::U_iterator u1 = qp1.u();
-    typename Quadratic_program2::U_iterator u2 = qp2.u();
-    for (int j=0; j<n; ++j, ++u1, ++u2)
-      if (*u1 != *u2) {
+      if ((*fu1 == true) && (*u1 != *u2)) {
 	std::cerr << "Equality test fails with u[" << j << "]: "
 		  << *u1 << " vs. " <<  *u2 << std::endl;
 	return_val = false;
       }
+    }
     // check d
     typename Quadratic_program1::D_iterator d1 = qp1.d();
     typename Quadratic_program2::D_iterator d2 = qp2.d();
@@ -299,19 +302,28 @@ namespace QP_functions_detail {
     typename P::C_iterator c = p.c();
     out << "COLUMNS\n";
     for (int j=0; j<n; ++j, ++c, ++a) {
-      if (!CGAL_NTS is_zero (*c))
+      // make sure that variable appears here even if it has only
+      // zero coefficients
+      bool written = false;
+      if (!CGAL_NTS is_zero (*c)) {
 	out << "  x" << j << "  obj  " << *c << "\n";
-      for (int i=0; i<m; ++i) { 
-	if (!CGAL_NTS is_zero ((*a)[i]))
-	  out << "  x" << j << "  c" << i << "  " << (*a)[i] << "\n";
+	written = true;
       }
+      for (int i=0; i<m; ++i) { 
+	if (!CGAL_NTS is_zero ((*a)[i])) {
+	  out << "  x" << j << "  c" << i << "  " << (*a)[i] << "\n";
+	  written = true;
+	}
+      }
+      if (!written)
+	out << "  x" << j << "  obj  " << 0 << "\n";
     }
  
     // RHS section:
     typename P::B_iterator b = p.b();
     out << "RHS\n";
     if (!CGAL_NTS is_zero (p.c0()))
-      out << "  rhs obj   -" << p.c0() << "\n";
+      out << "  rhs obj " << -p.c0() << "\n";
     for (int i=0; i<m; ++i, ++b)  
       if (!CGAL_NTS is_zero (*b))
 	out << "  rhs c" << i << "  " << *b << "\n";
diff --git a/QP_solver/include/CGAL/QP_solver/QP_solver_impl.h b/QP_solver/include/CGAL/QP_solver/QP_solver_impl.h
index 9e756b359e2..b48262eb603 100644
--- a/QP_solver/include/CGAL/QP_solver/QP_solver_impl.h
+++ b/QP_solver/include/CGAL/QP_solver/QP_solver_impl.h
@@ -161,8 +161,9 @@ solution_numerator( ) const
 	s +=  et2 * ET(qp_c[j]) * *j_it;
       z += d * s;
     }
-    // finally, add the constant term
-    return z += et2 * ET(qp_c0) * d * d;
+    // finally, add the constant term (phase II only)
+    if (is_phaseII) z += et2 * ET(qp_c0) * d * d;
+    return z;
 }
 
 // pivot step
@@ -929,7 +930,7 @@ test_mixed_bounds_dir_neg(int k, const ET& x_k, const ET& q_k,
             }
         }
     } else {                                        // check for upper bound
-        if ((q_k < et0) && (k < qp_n) && *(qp_fu+k)) {
+        if ((q_k > et0) && (k < qp_n) && *(qp_fu+k)) {
             ET  diff = (d * ET(qp_u[k])) - x_k;
             if (diff * q_min < d_min * q_k) {
                 i_min = k;
diff --git a/QP_solver/include/CGAL/QP_solver/QP_solver_validity_impl.h b/QP_solver/include/CGAL/QP_solver/QP_solver_validity_impl.h
index d76cdd11a90..59675389394 100644
--- a/QP_solver/include/CGAL/QP_solver/QP_solver_validity_impl.h
+++ b/QP_solver/include/CGAL/QP_solver/QP_solver_validity_impl.h
@@ -75,6 +75,7 @@ bool QP_solver<Q, ET, Tags>::is_valid() const
       {
 	const bool f = this->is_solution_feasible_for_auxiliary_problem();
 	const bool o = this->is_solution_optimal_for_auxiliary_problem();
+	const bool i = this->is_solution_infeasible();
 	const bool aux_positive = 
 	  ( (this->solution_numerator() > et0) && (d > et0) );
 	CGAL_qpe_debug {
@@ -85,9 +86,10 @@ bool QP_solver<Q, ET, Tags>::is_valid() const
 	  vout << "      is in phase I: " << is_phaseI << std::endl;
 	  vout << "       feasible_aux: " << f << std::endl;
 	  vout << "        optimal_aux: " << o << std::endl;
+	  vout << "        infeasible:  " << i << std::endl;
 	  vout << " obj. val. positive: " << aux_positive << std::endl;
 	}
-	return is_phaseI && f && o && aux_positive;
+	return is_phaseI && f && o && i && aux_positive;
       }
     case QP_UNBOUNDED:    
       {
@@ -127,16 +129,30 @@ bool QP_solver<Q, ET, Tags>::is_solution_feasible_for_auxiliary_problem() const
     if (*i_it < qp_n) {                 // original variable?
       if ((has_finite_lower_bound(*i_it) && (*v_it < lower_bound(*i_it) * d))||
 	  (has_finite_upper_bound(*i_it) && (*v_it > upper_bound(*i_it) * d)))
-        return false;
+	{
+	  CGAL_qpe_debug {
+	    vout4 << "variable " << *i_it << " out of bounds" << std::endl;
+	  }
+	  return false;
+	}
     } else                              // artificial variable?
       if (*v_it < et0)
-        return false;
+	{
+	  CGAL_qpe_debug {
+	    vout4 << "artificial variable negative" << std::endl;
+	  }
+	  return false;
+	}
   }
   
   // check nonegativity of slack variables (the basic ones suffice):
   for (Value_const_iterator v_it = x_B_S.begin(); v_it != x_B_S.end(); ++v_it)
-    if (*v_it < et0)
+    if (*v_it < et0) {
+      CGAL_qpe_debug {
+	vout4 << "slack variable out of bounds" << std::endl;
+      }
       return false;
+    }
   
   // Check whether the current solution is feasible for the auxiliary
   // problem.  For this, we use the fact that the auxiliary problem
@@ -192,11 +208,75 @@ bool QP_solver<Q, ET, Tags>::is_solution_feasible_for_auxiliary_problem() const
 
   // check equality (C1);
   for (int i=0; i<qp_m; ++i)
-    if (lhs_col[i] != ET(qp_b[i]) * d)
+    if (lhs_col[i] != ET(qp_b[i]) * d) {
+      CGAL_qpe_debug {
+	vout4 << "row " << i << " infeasible" << std::endl;
+      }
       return false;
+    }
   return true;
 }
 
+template < typename Q, typename ET, typename Tags >
+bool QP_solver<Q, ET, Tags>::is_solution_infeasible() const
+{
+  // checks whether we have a proof of infeasibilty according
+  // to Farkas Lemma. Namely, the system is infeasible iff
+  //     \tau^T = \lambda^T A  (A)
+  //     \lambda^T b - \sum{j: \tau_j < 0} \tau_j u_j  
+  //                 - \sum{j: \tau_j > 0} \tau_j l_j < 0  (B)
+  //      \lambda_i >= 0                for <=-constraints i,    (C)
+  //      \lambda_i <= 0                for >=-constraints i,    (D)
+ 
+  // get \lambda:
+  // Note: as the \lambda' we have access to is actual d times the \lambda
+  // we will not compute \tau but \tau * d.
+  Optimality_certificate_numerator_iterator 
+    itb = optimality_certificate_numerator_begin();
+  Optimality_certificate_numerator_iterator 
+    ite = optimality_certificate_numerator_end();
+  Optimality_certificate_iterator 
+    itq = optimality_certificate_begin();
+  Values lambda_prime;
+  int row = 0;
+  for (Optimality_certificate_numerator_iterator 
+	 it = itb; it != ite; ++it, ++itq, ++row) {
+    // simply check whether numerator/denominator gives correct value
+    CGAL_qpe_assertion (*it * (*itq).denominator() == (*itq).numerator() * d);
+    lambda_prime.push_back(*it);
+    // check sign (C) (D)
+    if (qp_r[row] == CGAL::SMALLER)
+      if (lambda_prime[row] < et0) return false;
+    if (qp_r[row] == CGAL::LARGER)
+      if (lambda_prime[row] > et0) return false;
+  }
+    
+  // compute \tau^T = \lambda^T A (A):
+  Values tau(qp_n,et0);
+  for (int col = 0; col < qp_n; ++col)
+    for (int i=0; i<qp_m; ++i)
+      tau[col] += lambda_prime[i] * ET((*(qp_A+col))[i]);
+
+  // check condition (B)
+  ET sum = et0;
+  // add \lambda^T b
+  for (int i=0; i<qp_m; ++i)
+    sum += lambda_prime[i] * ET(qp_b[i]);  // lambda carries factor of d
+  // subtract the remaining terms
+  for (int col = 0; col < qp_n; ++col) {
+    if (tau[col] < et0) {
+      CGAL_qpe_assertion(has_finite_upper_bound(col));
+      sum -= tau[col] * upper_bound(col); // tau carries factor of d
+    }
+    if (tau[col] > et0) {
+      CGAL_qpe_assertion(has_finite_lower_bound(col));
+      sum -= tau[col] * lower_bound(col); // tau carries factor of d
+    }
+  }
+  return sum < et0;
+}
+
+
 template < typename Q, typename ET, typename Tags >
 bool QP_solver<Q, ET, Tags>::is_value_correct() const
 {
@@ -443,8 +523,12 @@ bool QP_solver<Q, ET, Tags>::is_solution_feasible() const
     if (var != *it)
       return false; // original_variables_numerator_begin() inconsistent
     if ((has_finite_lower_bound(i) && var < lower_bound(i) * d) ||
-	(has_finite_upper_bound(i) && var > upper_bound(i) * d))
+	(has_finite_upper_bound(i) && var > upper_bound(i) * d)) {
+      CGAL_qpe_debug {
+	vout4 << "variable " << i << " out of bounds" << std::endl;
+      }
       return false;
+    }
 
     // compute A x times d:
     for (int j=0; j<qp_m; ++j)
@@ -456,8 +540,12 @@ bool QP_solver<Q, ET, Tags>::is_solution_feasible() const
     const ET rhs = ET(qp_b[row])*d;
     if ((qp_r[row] == CGAL::EQUAL         && lhs_col[row] != rhs) ||
 	(qp_r[row] == CGAL::SMALLER    && lhs_col[row] >  rhs) ||
-	(qp_r[row] == CGAL::LARGER && lhs_col[row] <  rhs))
+	(qp_r[row] == CGAL::LARGER && lhs_col[row] <  rhs)) {
+      CGAL_qpe_debug {
+	vout4 << "row " << row << " infeasible" << std::endl;
+      }
       return false;
+    }
   }
   
   return true;
diff --git a/QP_solver/test/QP_solver/master_mps_to_derivatives.cpp b/QP_solver/test/QP_solver/master_mps_to_derivatives.cpp
index 77ddd972ac2..7ff6314fb9f 100644
--- a/QP_solver/test/QP_solver/master_mps_to_derivatives.cpp
+++ b/QP_solver/test/QP_solver/master_mps_to_derivatives.cpp
@@ -191,7 +191,7 @@ void create_shifted_instance(CGAL::QP_from_mps
   for (int i=0; i<n; ++i) {
     for (int j=0; j<n; ++j)
       mvTD[i] += qp.d()[j][i] * v[j];
-    mvTD[i] *= -2;
+    mvTD[i] *= -1;
   }
 
   // output:
diff --git a/QP_solver/test/QP_solver/test_MPS.cpp b/QP_solver/test/QP_solver/test_MPS.cpp
index 61722f05555..7dbcc1a836e 100644
--- a/QP_solver/test/QP_solver/test_MPS.cpp
+++ b/QP_solver/test/QP_solver/test_MPS.cpp
@@ -25,7 +25,7 @@
 #ifndef CGAL_USE_GMP
 #include <CGAL/MP_Float.h>
 #else
-#include <CGAL/Gmpzf.h>
+#include <CGAL/Gmpz.h>
 #endif
 
 #include <CGAL/QP_solver.h>
@@ -49,11 +49,11 @@ int main(const int argNr,const char **args) {
   const int verbosity = argNr < 2? 1 : std::atoi(args[1]);
 
   // construct QP instance:
-  typedef double IT;
+  typedef int IT;
 #ifndef CGAL_USE_GMP
   typedef CGAL::MP_Float ET;
 #else
-  typedef CGAL::Gmpzf ET;
+  typedef CGAL::Gmpz ET;
 #endif
   typedef CGAL::QP_from_mps<IT, CGAL::Tag_false, CGAL::Tag_true> QP;
   QP qp(std::cin,true,verbosity);
@@ -69,7 +69,7 @@ int main(const int argNr,const char **args) {
     CGAL::print_quadratic_program (cout, qp);
     cout << std::endl;
   }
-  typedef CGAL::QP_solver_impl::QP_tags<CGAL::Tag_false,CGAL::Tag_true> Tags;
+  typedef CGAL::QP_solver_impl::QP_tags<CGAL::Tag_false,CGAL::Tag_false> Tags;
   CGAL::QP_pricing_strategy<QP, ET, Tags> *pricing_strategy =
     new CGAL::QP_full_exact_pricing<QP, ET, Tags>;
   typedef CGAL::QP_solver<QP, ET, Tags> Solver;