Merge pull request #11 from gergondet/topic/FixInfiniteLoop

src/CMakeLists.txt

@@ -63,7 +63,7 @@ endif()
 endif()
 
 if(${USE_F2C})
-  target_link_libraries(${PROJECT_NAME} libf2c.a)
+  target_link_libraries(${PROJECT_NAME} PUBLIC libf2c.a)
 endif()
 
 install(

src/QuadProg.cpp

@@ -32,7 +32,9 @@ QuadProgCommon::QuadProgCommon():
 	X_(),
 	fail_(0),
 	iact_(),
-	iter_(2)
+	iter_(2),
+	tol_(0.0),
+	maxiter_(0)
 {
 }
 
@@ -48,6 +50,34 @@ int QuadProgCommon::fail() const
 	return fail_;
 }
 
+int QuadProgCommon::maxiter() const
+{
+	return maxiter_;
+}
+
+void QuadProgCommon::maxiter(int maxiter)
+{
+	if(maxiter < 0)
+	{
+		throw std::domain_error("Maximum iteration count must be >= 0");
+	}
+	maxiter_ = maxiter;
+}
+
+double QuadProgCommon::tolerance() const
+{
+	return tol_;
+}
+
+void QuadProgCommon::tolerance(double tol)
+{
+	if(tol < 0.0)
+	{
+		tolerance(-tol);
+		return;
+	}
+	tol_ = tol;
+}
 
 const VectorXd& QuadProgCommon::result() const
 {
@@ -146,7 +176,7 @@ bool QuadProgDense::solve(const Ref<const MatrixXd>& Q, const Ref<const VectorXd
 
 	qpgen2_(Q_.data(), C_.data(), &fddmat, &n, X_.data(), &crval,
 		A_.data(), B_.data(), &fdamat, &q, &meq, iact_.data(), &nact,
-		iter_.data(), work_.data(), &fail_);
+		iter_.data(), work_.data(), &fail_, &tol_, &maxiter_);
 
 	return fail_ == 0;
 }
@@ -226,7 +256,7 @@ bool QuadProgSparse::solve(const Ref<const MatrixXd>& Q, const Ref<const VectorX
 
 	qpgen1_(Q_.data(), C_.data(), &fddmat, &n, X_.data(), &crval,
 		A_.data(), iA_.data(), B_.data(), &fdamat, &q, &meq, iact_.data(), &nact,
-		iter_.data(), work_.data(), &fail_);
+		iter_.data(), work_.data(), &fail_, &tol_, &maxiter_);
 
 	return fail_ == 0;
 }

src/QuadProg.h

@@ -29,12 +29,12 @@ namespace Eigen
 extern "C" int qpgen1_(double* dmat, double* dvec, const int* fddmat,
 	const int* n, double* sol, double* crval, double* amat, const int* iamat,
 	double* bvec, const int* fdamat, const int* q, const int* meq, int* iact,
-	int* nact, int* iter, double* work, const int* ierr);
+	int* nact, int* iter, double* work, const int* ierr, double* tol, int* maxiter);
 
 extern "C" int qpgen2_(double* dmat, double* dvec, const int* fddmat,
 	const int* n, double* sol, double* crval, double* amat, double* bvec,
 	const int* fdamat, const int* q, const int* meq, int* iact, int* nact,
-	int* iter, double* work, const int* ierr);
+	int* iter, double* work, const int* ierr, double* tol, int* maxiter);
 
 /** Common method for Quadprog solver classes.
  *
@@ -71,6 +71,32 @@ class QuadProgCommon
      */
 	EIGEN_QUADPROG_API const VectorXd& result() const;
 
+    /** Maximum iteration count
+     *
+     * Defaults to max(50, 5 * (nrvar + nreq + nrineq) if 0
+     *
+     */
+	EIGEN_QUADPROG_API int maxiter() const;
+
+    /** Set the maximum iteration count
+     *
+     */
+	EIGEN_QUADPROG_API void maxiter(int maxiter);
+
+    /** Constraint violation tolerance used by the solver
+     *
+     * Throw if maxiter < 0
+     *
+     */
+	EIGEN_QUADPROG_API double tolerance() const;
+
+    /** Set the constraint violation tolerance used by the solver
+     *
+     * If tol < 0, act as tolerance(-tol)
+     *
+     */
+	EIGEN_QUADPROG_API void tolerance(double tol);
+
     /** Set problem dimensions.
      *
      * \param nrvar Dimension \f$n\f$ of optimization vector \f$x \in
@@ -103,6 +129,8 @@ class QuadProgCommon
                       deleted after they became active */
     VectorXd work_; /**< Working space vector with length at least
                       \f$2n+r(r+5)/2+2q+1\f$ where \f$r=\min(n,q)\f$ */
+	double tol_; /**< Constraint violation tolerance */
+	int maxiter_; /**< Maximum iteration count */
 };
 
 /** Dense quadratic program.

src/QuadProg/c/solve.QP.c

@@ -1,4 +1,4 @@
-/* f/solve.QP.f -- translated by f2c (version 20100827).
+/* solve.QP.f -- translated by f2c (version 20160102).
    You must link the resulting object file with libf2c:
 	on Microsoft Windows system, link with libf2c.lib;
 	on Linux or Unix systems, link with .../path/to/libf2c.a -lm
@@ -72,6 +72,8 @@
 /*            ierr =  0, we have to decompose D */
 /*            ierr != 0, D is already decomposed into D=R^TR and we were */
 /*                       given R^{-1}. */
+/*  tol    scalar, constraint violation tolerance allowed by the solver */
+/*  maxiter integer, maximum iteration count */
 
 /*  Output parameter: */
 /*  sol   nx1 the final solution (x in the notation above) */
@@ -87,6 +89,7 @@
 /*           ierr = 1, the minimization problem has no solution */
 /*           ierr = 2, problems with decomposing D, in this case sol */
 /*                     contains garbage!! */
+/*           ierr = 3, max iterations reached */
 
 /*  Working space: */
 /*  work  vector with length at least 2*n+r*(r+5)/2 + 2*q +1 */
@@ -96,7 +99,7 @@
 	fddmat, integer *n, doublereal *sol, doublereal *crval, doublereal *
 	amat, doublereal *bvec, integer *fdamat, integer *q, integer *meq, 
 	integer *iact, integer *nact, integer *iter, doublereal *work, 
-	integer *ierr)
+	integer *ierr, doublereal *tol, integer *maxiter)
 {
     /* System generated locals */
     integer dmat_dim1, dmat_offset, amat_dim1, amat_offset, i__1, i__2;
@@ -136,6 +139,11 @@
     /* Function Body */
     r__ = min(*n,*q);
     l = (*n << 1) + r__ * (r__ + 5) / 2 + (*q << 1) + 1;
+    if (*maxiter == 0) {
+/* Computing MAX */
+	i__1 = 50, i__2 = (*n + *q) * 5;
+	*maxiter = max(i__1,i__2);
+    }
 
 /* store the initial dvec to calculate below the unconstrained minima of */
 /* the critical value. */
@@ -244,6 +252,10 @@
 /* start a new iteration */
 
     ++iter[1];
+    if (iter[1] > *maxiter) {
+	*ierr = 3;
+	goto L999;
+    }
 
 /* calculate all constraints and check which are still violated */
 /* for the equality constraints we have to check whether the normal */
@@ -291,7 +303,7 @@
 /* take always the first constraint which is violated. ;-) */
 
     nvl = 0;
-    temp = 0.;
+    temp = -(*tol);
     i__1 = *q;
     for (i__ = 1; i__ <= i__1; ++i__) {
 	if (work[iwsv + i__] < temp * work[iwnbv + i__]) {

src/QuadProg/c/solve.QP.compact.c

@@ -1,4 +1,4 @@
-/* f/solve.QP.compact.f -- translated by f2c (version 20100827).
+/* solve.QP.compact.f -- translated by f2c (version 20160102).
    You must link the resulting object file with libf2c:
 	on Microsoft Windows system, link with libf2c.lib;
 	on Linux or Unix systems, link with .../path/to/libf2c.a -lm
@@ -81,6 +81,8 @@
 /*            ierr =  0, we have to decompose D */
 /*            ierr != 0, D is already decomposed into D=R^TR and we were */
 /*                       given R^{-1}. */
+/*  tol    scalar, constraint violation tolerance allowed by the solver */
+/*  maxiter  integer, maximum iteration count */
 
 /*  Output parameter: */
 /*  sol   nx1 the final solution (x in the notation above) */
@@ -96,6 +98,7 @@
 /*           ierr = 1, the minimization problem has no solution */
 /*           ierr = 2, problems with decomposing D, in this case sol */
 /*                     contains garbage!! */
+/*           ierr = 3, max iterations reached */
 
 /*  Working space: */
 /*  work  vector with length at least 2*n+r*(r+5)/2 + 2*q +1 */
@@ -105,7 +108,7 @@
 	fddmat, integer *n, doublereal *sol, doublereal *crval, doublereal *
 	amat, integer *iamat, doublereal *bvec, integer *fdamat, integer *q, 
 	integer *meq, integer *iact, integer *nact, integer *iter, doublereal 
-	*work, integer *ierr)
+	*work, integer *ierr, doublereal *tol, integer *maxiter)
 {
     /* System generated locals */
     integer iamat_dim1, iamat_offset, dmat_dim1, dmat_offset, amat_dim1, 
@@ -149,6 +152,11 @@
     /* Function Body */
     r__ = min(*n,*q);
     l = (*n << 1) + r__ * (r__ + 5) / 2 + (*q << 1) + 1;
+    if (*maxiter == 0) {
+/* Computing MAX */
+	i__1 = 50, i__2 = (*n + *q) * 5;
+	*maxiter = max(i__1,i__2);
+    }
 
 /* store the initial dvec to calculate below the unconstrained minima of */
 /* the critical value. */
@@ -257,6 +265,10 @@
 /* start a new iteration */
 
     ++iter[1];
+    if (iter[1] > *maxiter) {
+	*ierr = 3;
+	goto L999;
+    }
 
 /* calculate all constraints and check which are still violated */
 /* for the equality constraints we have to check whether the normal */
@@ -305,7 +317,7 @@
 /* take always the first constraint which is violated. ;-) */
 
     nvl = 0;
-    temp = 0.;
+    temp = -(*tol);
     i__1 = *q;
     for (i__ = 1; i__ <= i__1; ++i__) {
 	if (work[iwsv + i__] < temp * work[iwnbv + i__]) {

src/QuadProg/f/solve.QP.compact.f

@@ -67,6 +67,8 @@
 c            ierr =  0, we have to decompose D
 c            ierr != 0, D is already decomposed into D=R^TR and we were
 c                       given R^{-1}.
+c  tol    scalar, constraint violation tolerance allowed by the solver
+c  maxiter  integer, maximum iteration count
 c
 c  Output parameter:
 c  sol   nx1 the final solution (x in the notation above)
@@ -82,24 +84,29 @@
 c           ierr = 1, the minimization problem has no solution
 c           ierr = 2, problems with decomposing D, in this case sol
 c                     contains garbage!!
+c           ierr = 3, max iterations reached
 c
 c  Working space:
 c  work  vector with length at least 2*n+r*(r+5)/2 + 2*q +1
 c        where r=min(n,q)
 c
       subroutine qpgen1(dmat, dvec, fddmat, n, sol, crval, amat, iamat,
-     *     bvec, fdamat, q, meq, iact, nact, iter, work, ierr)  
+     *     bvec, fdamat, q, meq, iact, nact, iter, work, ierr, tol,
+     *     maxiter)
       implicit none
       integer n, i, j, l, l1, fdamat, fddmat,
      *     info, q, iamat(fdamat+1,*), iact(*), iter(*), it1,
      *     ierr, nact, iwzv, iwrv, iwrm, iwsv, iwuv, nvl,
-     *     r, iwnbv, meq
+     *     r, iwnbv, meq, maxiter
       double precision dmat(fddmat,*), dvec(*),sol(*), bvec(*),
      *     work(*), temp, sum, t1, tt, gc, gs, crval,
-     *     nu, amat(fdamat,*)
+     *     nu, amat(fdamat,*), tol
       logical t1inf, t2min
       r = min(n,q)
       l = 2*n + (r*(r+5))/2 + 2*q + 1
+      if( maxiter .EQ. 0) then
+        maxiter = max(50, 5 * (n + q))
+      endif
 c 
 c store the initial dvec to calculate below the unconstrained minima of
 c the critical value.
@@ -186,6 +193,10 @@ subroutine qpgen1(dmat, dvec, fddmat, n, sol, crval, amat, iamat,
 c start a new iteration      
 c
       iter(1) = iter(1)+1
+      if (iter(1) .GT. maxiter) then
+        ierr = 3
+        goto 999
+      endif
 c
 c calculate all constraints and check which are still violated
 c for the equality constraints we have to check whether the normal
@@ -225,7 +236,7 @@ subroutine qpgen1(dmat, dvec, fddmat, n, sol, crval, amat, iamat,
 c take always the first constraint which is violated. ;-)
 c
       nvl = 0 
-      temp = 0.d0
+      temp = -tol
       do 71 i=1,q
          if (work(iwsv+i) .LT. temp*work(iwnbv+i)) then
             nvl = i

src/QuadProg/f/solve.QP.f

@@ -58,6 +58,8 @@
 c            ierr =  0, we have to decompose D
 c            ierr != 0, D is already decomposed into D=R^TR and we were
 c                       given R^{-1}.
+c  tol    scalar, constraint violation tolerance allowed by the solver
+c  maxiter integer, maximum iteration count
 c
 c  Output parameter:
 c  sol   nx1 the final solution (x in the notation above)
@@ -73,24 +75,29 @@
 c           ierr = 1, the minimization problem has no solution
 c           ierr = 2, problems with decomposing D, in this case sol
 c                     contains garbage!!
+c           ierr = 3, max iterations reached
 c
 c  Working space:
 c  work  vector with length at least 2*n+r*(r+5)/2 + 2*q +1
 c        where r=min(n,q)
 c
-      subroutine qpgen2(dmat, dvec, fddmat, n, sol, crval, amat, 
-     *     bvec, fdamat, q, meq, iact, nact, iter, work, ierr)  
+      subroutine qpgen2(dmat, dvec, fddmat, n, sol, crval, amat,
+     *     bvec, fdamat, q, meq, iact, nact, iter, work, ierr, tol,
+     *     maxiter)
       implicit none
       integer n, i, j, l, l1, fdamat, fddmat,
      *     info, q, iact(*), iter(*), it1,
      *     ierr, nact, iwzv, iwrv, iwrm, iwsv, iwuv, nvl,
-     *     r, iwnbv, meq
+     *     r, iwnbv, meq, maxiter
       double precision dmat(fddmat,*), dvec(*),sol(*), bvec(*),
      *     work(*), temp, sum, t1, tt, gc, gs, crval,
-     *     nu, amat(fdamat,*)
+     *     nu, amat(fdamat,*), tol
       logical t1inf, t2min
       r = min(n,q)
       l = 2*n + (r*(r+5))/2 + 2*q + 1
+      if( maxiter .EQ. 0) then
+        maxiter = max(50, 5 * (n + q))
+      endif
 c 
 c store the initial dvec to calculate below the unconstrained minima of
 c the critical value.
@@ -177,6 +184,10 @@ subroutine qpgen2(dmat, dvec, fddmat, n, sol, crval, amat,
 c start a new iteration      
 c
       iter(1) = iter(1)+1
+      if (iter(1) .GT. maxiter) then
+        ierr = 3
+        goto 999
+      endif
 c
 c calculate all constraints and check which are still violated
 c for the equality constraints we have to check whether the normal
@@ -216,7 +227,7 @@ subroutine qpgen2(dmat, dvec, fddmat, n, sol, crval, amat,
 c take always the first constraint which is violated. ;-)
 c
       nvl = 0 
-      temp = 0.d0
+      temp = -tol
       do 71 i=1,q
          if (work(iwsv+i) .LT. temp*work(iwnbv+i)) then
             nvl = i