From: stevenj <stevenj@alum.mit.edu>
Date: Tue, 19 Aug 2008 00:05:30 +0000 (-0400)
Subject: more accurate dual minimization for small u
X-Git-Url: http://www.chiark.greenend.org.uk/ucgi/~ianmdlvl/git?a=commitdiff_plain;h=ab8009c70a44e3debb30911ccad09f9cf5731d06;p=nlopt.git

more accurate dual minimization for small u

darcs-hash:20080819000530-c8de0-f967384aa86b75bdac46769c2929ba51a1d4bca1.gz
---

diff --git a/mma/mma.c b/mma/mma.c
index f1b72d1..7618306 100644
--- a/mma/mma.c
+++ b/mma/mma.c
@@ -74,7 +74,12 @@ static double dual_func(int m, const double *y, double *grad, void *d_)
 	       v += (fabs(dfcdx[i*n + j]) * sigma[j] + 0.5 * rhoc[i]) * y[i];
 	  }
 	  u *= (sigma2 = sqr(sigma[j]));
-	  dx = u==0 ? 0 : (v/u)*sigma2 * (-1 + sqrt(1 - sqr(u/(v*sigma[j]))));
+	  if (fabs(u) < 1e-3 * (v*sigma[j])) { /* Taylor exp. for small u */
+	       double a = u / (v*sigma[j]);
+	       dx = -sigma[j] * (0.5 * a + 0.125 * a*a*a);
+	  }
+	  else
+	       dx = (v/u)*sigma2 * (-1 + sqrt(1 - sqr(u/(v*sigma[j]))));
 	  xcur[j] = x[j] + dx;
 	  if (xcur[j] > ub[j]) xcur[j] = ub[j];
 	  else if (xcur[j] < lb[j]) xcur[j] = lb[j];