Add the option to set the maximum number of iterations

src/QuadProg.cpp

@@ -33,7 +33,8 @@ QuadProgCommon::QuadProgCommon():
 	fail_(0),
 	iact_(),
 	iter_(2),
-	tol_(0.0)
+	tol_(0.0),
+	maxiter_(0)
 {
 }
 
@@ -49,6 +50,20 @@ 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_;
@@ -161,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_, &tol_);
+		iter_.data(), work_.data(), &fail_, &tol_, &maxiter_);
 
 	return fail_ == 0;
 }
@@ -241,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_, &tol_);
+		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, double* tol);
+	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, double* tol);
+	int* iter, double* work, const int* ierr, double* tol, int* maxiter);
 
 /** Common method for Quadprog solver classes.
  *
@@ -71,7 +71,21 @@ 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;
@@ -116,6 +130,7 @@ class QuadProgCommon
     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

@@ -73,6 +73,7 @@
 /*            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) */
@@ -98,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, doublereal *tol)
+	integer *ierr, doublereal *tol, integer *maxiter)
 {
     /* System generated locals */
     integer dmat_dim1, dmat_offset, amat_dim1, amat_offset, i__1, i__2;
@@ -119,7 +120,7 @@
     extern /* Subroutine */ int dpofa_(doublereal *, integer *, integer *, 
 	    integer *), dpori_(doublereal *, integer *, integer *), dposl_(
 	    doublereal *, integer *, integer *, doublereal *);
-    static integer iwnbv, maxiter;
+    static integer iwnbv;
 
     /* Parameter adjustments */
     --dvec;
@@ -138,9 +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);
+	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. */
@@ -249,7 +252,7 @@
 /* start a new iteration */
 
     ++iter[1];
-    if (iter[1] > maxiter) {
+    if (iter[1] > *maxiter) {
 	*ierr = 3;
 	goto L999;
     }

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

@@ -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) */
@@ -97,7 +99,6 @@
 /*           ierr = 2, problems with decomposing D, in this case sol */
 /*                     contains garbage!! */
 /*           ierr = 3, max iterations reached */
-/*  tol    scalar, constraint violation tolerance allowed by the solver */
 
 /*  Working space: */
 /*  work  vector with length at least 2*n+r*(r+5)/2 + 2*q +1 */
@@ -107,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, doublereal *tol)
+	*work, integer *ierr, doublereal *tol, integer *maxiter)
 {
     /* System generated locals */
     integer iamat_dim1, iamat_offset, dmat_dim1, dmat_offset, amat_dim1, 
@@ -129,7 +130,7 @@
     extern /* Subroutine */ int dpofa_(doublereal *, integer *, integer *, 
 	    integer *), dpori_(doublereal *, integer *, integer *), dposl_(
 	    doublereal *, integer *, integer *, doublereal *);
-    static integer iwnbv, maxiter;
+    static integer iwnbv;
 
     /* Parameter adjustments */
     --dvec;
@@ -151,9 +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);
+	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. */
@@ -262,7 +265,7 @@
 /* start a new iteration */
 
     ++iter[1];
-    if (iter[1] > maxiter) {
+    if (iter[1] > *maxiter) {
 	*ierr = 3;
 	goto L999;
     }

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)
@@ -83,14 +85,14 @@
 c           ierr = 2, problems with decomposing D, in this case sol
 c                     contains garbage!!
 c           ierr = 3, max iterations reached
-c  tol    scalar, constraint violation tolerance allowed by the solver
 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, tol)
+     *     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,
@@ -102,7 +104,9 @@ subroutine qpgen1(dmat, dvec, fddmat, n, sol, crval, amat, iamat,
       logical t1inf, t2min
       r = min(n,q)
       l = 2*n + (r*(r+5))/2 + 2*q + 1
-      maxiter = max(50, 5 * (n + q))
+      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.

src/QuadProg/f/solve.QP.f

@@ -59,6 +59,7 @@
 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)
@@ -81,7 +82,8 @@
 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, tol)
+     *     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,
@@ -93,7 +95,9 @@ subroutine qpgen2(dmat, dvec, fddmat, n, sol, crval, amat,
       logical t1inf, t2min
       r = min(n,q)
       l = 2*n + (r*(r+5))/2 + 2*q + 1
-      maxiter = max(50, 5 * (n + q))
+      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.