bug fixes to slsqp, now seems to work

darcs-hash:20100712223028-c8de0-5102cafe2ad38f4e100fb986cb3085cb0e7f2aba.gz

diff --git a/slsqp/README b/slsqp/README
index 4dd7593c8018a644745abbc2bf9d29e1f221c7fd..a24b0339ce48681c7dd1e5a8d7a1a50f5e23b740 100644 (file)
--- a/slsqp/README
+++ b/slsqp/README
@@ -21,11 +21,16 @@ cleaned up.  It was modified to be re-entrant, preserving the
 reverse-communication interface but explicitly saving the state in a
 data structure.  The reverse-communication interface was wrapped with
 an NLopt-style inteface, with NLopt stopping conditions.  The inexact
-line search was modified to evaluate the functions including
-gradients, since this removes the need to evaluate the
-function+gradient a second time for the same point when the inexact
-line search concludes (hopefully after one or two steps), since
+line search was modified to evaluate the functions including gradients
+for the first step, since this removes the need to evaluate the
+function+gradient a second time for the same point in the common case
+when the inexact line search concludes after a single step, since
 NLopt's interface combines the function and gradient computations.
+Since roundoff errors sometimes pushed SLSQP's parameters slightly
+outside the bound constraints (not allowed by NLopt), we added checks
+to force the parameters within the bounds.  Fixed a bug in LSEI (use
+of uninitialized variables) for the case where the number of equality
+constraints equals the dimension of the problem.
 
 The exact line-search option is currently disabled; if we want to
 re-enable this (although exact line-search is usually overkill in

diff --git a/slsqp/slsqp.c b/slsqp/slsqp.c
index 832ac80d08937fda839e8bf5b4e6f7d589f7e6b0..632e2615202eb234d9ae7651e105a4e1abb8733e 100644 (file)
--- a/slsqp/slsqp.c
+++ b/slsqp/slsqp.c
@@ -1145,7 +1145,7 @@ static void lsei_(double *c__, double *d__, double *e,
     }
     *mode = 1;
     w[mc1] = 0.0;
-    i__2 = *mg - *mc;
+    i__2 = *mg; // BUGFIX for *mc == *n: changed from *mg - *mc, SGJ 2010
     dcopy___(&i__2, &w[mc1], 0, &w[mc1], 1);
     if (*mc == *n) {
        goto L50;
@@ -1413,6 +1413,15 @@ static void lsq_(int *m, int *meq, int *n, int *nl,
        dcopy___(m, &w[iw], 1, &y[1], 1);
        dcopy___(&n3, &w[iw + *m], 1, &y[*m + 1], 1);
        dcopy___(&n3, &w[iw + *m + *n], 1, &y[*m + n3 + 1], 1);
+
+       /* SGJ, 2010: make sure bound constraints are satisfied, since
+          roundoff error sometimes causes slight violations and
+          NLopt guarantees that bounds are strictly obeyed */
+       i__1 = *n;
+       for (i__ = 1; i__ <= i__1; ++i__) {
+            if (x[i__] < xl[i__]) x[i__] = xl[i__];
+            else if (x[i__] > xu[i__]) x[i__] = xu[i__];
+       }
     }
 /*   END OF SUBROUTINE LSQ */
 } /* lsq_ */
@@ -1901,6 +1910,7 @@ L130:
     h4 = one;
     lsq_(m, meq, n, &n3, la, &l[1], &g[1], &a[a_offset], &c__[1], &u[1], &v[1]
            , &s[1], &r__[1], &w[1], &iw[1], mode);
+
 /*   AUGMENTED PROBLEM FOR INCONSISTENT LINEARIZATION */
     if (*mode == 6) {
        if (*n == *meq) {
@@ -2010,13 +2020,19 @@ L190:
     dscal_sl__(n, &alpha, &s[1], 1);
     dcopy___(n, &x0[1], 1, &x[1], 1);
     daxpy_sl__(n, &one, &s[1], 1, &x[1], 1);
-    /* SGJ change, 2010: since we should expect the inexact linesearch to often succeed on the
-       first try or two, there is no point in not computing the gradient information (especially
-       since gradients are cheap).  Hence we pass a special mode -2, that tells the caller
-       to evaluate the function and gradient values (like mode -1), but also tells the caller
-       that the subsequent mode -1 call can return without evaluating the function (since
-       the subsequent mode -1 call will always be for the same point as the last mode -2 call) */
-    *mode = -2;
+    
+    /* SGJ 2010: ensure roundoff doesn't push us past bound constraints */
+    i__1 = *n; for (i__ = 1; i__ <= i__1; ++i__) {
+        if (x[i__] < xl[i__]) x[i__] = xl[i__];
+        else if (x[i__] > xu[i__]) x[i__] = xu[i__];
+    }
+
+    /* SGJ 2010: optimizing for the common case where the inexact line
+       search succeeds in one step, use special mode = -2 here to
+       eliminate a a subsequent unnecessary mode = -1 call, at the 
+       expense of extra gradient evaluations when more than one inexact
+       line-search step is required */
+    *mode = line == 1 ? -2 : 1;
     goto L330;
 L200:
     if (h1 <= h3 / ten || line > 10) {
@@ -2414,7 +2430,7 @@ nlopt_result nlopt_slsqp(unsigned n, nlopt_func f, void *f_data,
 {
      slsqpb_state state = {0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0};
      unsigned mtot = nlopt_count_constraints(m, fc);
-     unsigned ptot = nlopt_count_constraints(p, fc);
+     unsigned ptot = nlopt_count_constraints(p, h);
      double *work, *cgrad, *c, *grad, *w, 
          fcur, *xcur, fprev, *xprev, *cgradtmp;
      int mpi = (int) (mtot + ptot), pi = (int) ptot,  ni = (int) n, mpi1 = mpi > 0 ? mpi : 1;
@@ -2424,8 +2440,8 @@ nlopt_result nlopt_slsqp(unsigned n, nlopt_func f, void *f_data,
      int iter = 0; /* tell sqsqp to ignore this check, since we check evaluation counts ourselves */
      unsigned i, ii;
      nlopt_result ret = NLOPT_SUCCESS;
-     int feasible = 0, feasible_cur;
-     double infeasibility = HUGE_VAL, infeasibility_cur;
+     int feasible, feasible_cur;
+     double infeasibility = HUGE_VAL, infeasibility_cur = HUGE_VAL;
      unsigned max_cdim;
      
      max_cdim = max(nlopt_max_constraint_dim(m, fc),
@@ -2434,15 +2450,15 @@ nlopt_result nlopt_slsqp(unsigned n, nlopt_func f, void *f_data,
 
 #define U(n) ((unsigned) (n))
      work = (double *) malloc(sizeof(double) * (U(mpi1) * (n + 1) 
-                                               + U(mtot+ptot) 
-                                               + n + n + n + max_cdim*n
+                                               + U(mpi) 
+                                               + n+1 + n + n + max_cdim*n
                                                + U(len_w))
                              + sizeof(int) * U(len_jw));
      if (!work) return NLOPT_OUT_OF_MEMORY;
      cgrad = work;
      c = cgrad + U(mpi1) * (n + 1);
-     grad = c + (mtot+ptot);
-     xcur = grad + n;
+     grad = c + mpi;
+     xcur = grad + n+1;
      xprev = xcur + n;
      cgradtmp = xprev + n;
      w = cgradtmp + max_cdim*n;
@@ -2451,7 +2467,9 @@ nlopt_result nlopt_slsqp(unsigned n, nlopt_func f, void *f_data,
      memcpy(xcur, x, sizeof(double) * n);
      memcpy(xprev, x, sizeof(double) * n);
      fprev = fcur = *minf = HUGE_VAL;
-     feasible = 0;
+     feasible = feasible_cur = 0;
+
+     goto eval_f_and_grad; /* eval before calling slsqp the first time */
 
      do {
          slsqp(&mpi, &pi, &mpi1, &ni,
@@ -2465,6 +2483,7 @@ nlopt_result nlopt_slsqp(unsigned n, nlopt_func f, void *f_data,
              case -1:  /* objective & gradient evaluation */
                   if (prev_mode == -2) break; /* just evaluated this point */
              case -2:
+                  eval_f_and_grad:
                   feasible_cur = 1; infeasibility_cur = 0;
                   fcur = f(n, xcur, grad, f_data);
                   if (nlopt_stop_forced(stop)) { 
@@ -2537,7 +2556,7 @@ nlopt_result nlopt_slsqp(unsigned n, nlopt_func f, void *f_data,
              case 8: /* positive directional derivative for linesearch */
                   /* relaxed convergence check for a feasible_cur point,
                      as in the SLSQP code */
-                  ret = NLOPT_FAILURE;
+                  ret = NLOPT_ROUNDOFF_LIMITED; /* usually why deriv>0 */
                   if (feasible_cur) {
                        double save_ftol_rel = stop->ftol_rel;
                        double save_ftol_abs = stop->ftol_abs;
@@ -2570,6 +2589,7 @@ nlopt_result nlopt_slsqp(unsigned n, nlopt_func f, void *f_data,
          if ((fcur < *minf && (feasible_cur || !feasible))
              || (!feasible && infeasibility_cur < infeasibility)) {
               *minf = fcur;
+              feasible = feasible_cur;
               infeasibility = infeasibility_cur;
               memcpy(x, xcur, sizeof(double)*n);
          }