more accurate dual minimization for small u

darcs-hash:20080819000530-c8de0-f967384aa86b75bdac46769c2929ba51a1d4bca1.gz

diff --git a/mma/mma.c b/mma/mma.c
index f1b72d167be206ba57f47731b8c77cab4d190ba3..76183061a6d519b242169622ea60d4f445f0374b 100644 (file)
--- 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];