chiark
/
gitweb
/
~ianmdlvl
/
nlopt.git
/ commitdiff
commit
grep
author
committer
pickaxe
?
search:
re
summary
|
shortlog
|
log
|
commit
| commitdiff |
tree
raw
|
patch
|
inline
| side by side (parent:
5131a48
)
use reciprocal quadratic formula, which is accurate everywhere, rather than switching...
author
stevenj
<stevenj@alum.mit.edu>
Sat, 26 Nov 2011 19:45:06 +0000
(14:45 -0500)
committer
stevenj
<stevenj@alum.mit.edu>
Sat, 26 Nov 2011 19:45:06 +0000
(14:45 -0500)
Ignore-this:
61a4e0cf96604cc766bad8540402b95f
darcs-hash:
20111126194506
-c8de0-
3099879e86db6a94cba3ebd07aad2ccbfc968ff3
.gz
mma/mma.c
patch
|
blob
|
history
diff --git
a/mma/mma.c
b/mma/mma.c
index 03638152d4ff00bd64b64ce1a9e6572fe9aae177..249b892e0224aebfaecd55e6cbfb4018fc715877 100644
(file)
--- 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))
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];
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]));
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];
xcur[j] = x[j] + dx;
if (xcur[j] > ub[j]) xcur[j] = ub[j];
else if (xcur[j] < lb[j]) xcur[j] = lb[j];