From a96849e95ad3c5892e3357f4823a50ff83112891 Mon Sep 17 00:00:00 2001 From: stevenj Date: Sat, 26 Nov 2011 14:45:06 -0500 Subject: [PATCH] use reciprocal quadratic formula, which is accurate everywhere, rather than switching between quadratic formula and Taylor expansion Ignore-this: 61a4e0cf96604cc766bad8540402b95f darcs-hash:20111126194506-c8de0-3099879e86db6a94cba3ebd07aad2ccbfc968ff3.gz --- mma/mma.c | 11 ++++------- 1 file changed, 4 insertions(+), 7 deletions(-) diff --git a/mma/mma.c b/mma/mma.c index 0363815..249b892 100644 --- a/mma/mma.c +++ b/mma/mma.c @@ -94,7 +94,9 @@ static double dual_func(unsigned m, const double *y, double *grad, void *d_) and it follows that the only dx solution with |dx| <= sigma is given by: (v/u) sigma^2 (-1 + sqrt(1 - (u / v sigma)^2)) - (which goes to zero as u -> 0). */ + = (u/v) / (-1 - sqrt(1 - (u / v sigma)^2)) + (which goes to zero as u -> 0). The latter expression + is less susceptible to roundoff error. */ if (sigma[j] == 0) { /* special case for lb[i] == ub[i] dims, dx=0 */ xcur[j] = x[j]; @@ -108,12 +110,7 @@ static double dual_func(unsigned 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])); - 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(fabs(1 - sqr(u/(v*sigma[j]))))); + dx = (u/v) / (-1 - sqrt(fabs(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]; -- 2.30.2