Mercurial > octave

changeset 21088:323e92c4589f stable

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
fzero.m: Correctly choose tolerance (eps) based on class of fun and X0 (bug #46658). * fzero.m: Check whether X0 or the result of the evaluation of the function fun results in a variable of class "single". Use an eps which matches X0 or the function evaluation.
author Rik <rik@octave.org>
date Sat, 16 Jan 2016 18:14:35 -0800
parents 8fdf1f991a32
children 2d7f658daa58 0cfcc178432b
files scripts/optimization/fzero.m
diffstat 1 files changed, 10 insertions(+), 5 deletions(-) [+]
line wrap: on
line diff
--- a/scripts/optimization/fzero.m	Tue Jan 05 12:05:01 2016 -0800
+++ b/scripts/optimization/fzero.m	Sat Jan 16 18:14:35 2016 -0800
@@ -137,7 +137,6 @@
   nfev = 0;
 
   x = fval = a = fa = b = fb = NaN;
-  eps = eps (class (x0));
 
   ## Prepare...
   a = x0(1);
@@ -194,6 +193,12 @@
     u = b; fu = fb;
   endif
 
+  if (isa (x0, "single") || isa (fa, "single"))
+    macheps = eps ("single");
+  else
+    macheps = eps ("double");
+  endif
+
   d = e = u;
   fd = fe = fu;
   mba = mu*(b - a);
@@ -201,7 +206,7 @@
     switch (itype)
       case 1
         ## The initial test.
-        if (b - a <= 2*(2 * abs (u) * eps + tolx))
+        if (b - a <= 2*(2 * abs (u) * macheps + tolx))
           x = u; fval = fu;
           info = 1;
           break;
@@ -266,7 +271,7 @@
     endswitch
 
     ## Don't let c come too close to a or b.
-    delta = 2*0.7*(2 * abs (u) * eps + tolx);
+    delta = 2*0.7*(2 * abs (u) * macheps + tolx);
     if ((b - a) <= 2*delta)
       c = (a + b)/2;
     else
@@ -320,7 +325,7 @@
     else
       u = b; fu = fb;
     endif
-    if (b - a <= 2*(2 * abs (u) * eps + tolx))
+    if (b - a <= 2*(2 * abs (u) * macheps + tolx))
       info = 1;
       break;
     endif
@@ -337,7 +342,7 @@
   ## Check solution for a singularity by examining slope
   if (info == 1)
     if ((b - a) != 0
-        && abs ((fb - fa)/(b - a) / slope0) > max (1e6, 0.5/(eps+tolx)))
+        && abs ((fb - fa)/(b - a) / slope0) > max (1e6, 0.5/(macheps+tolx)))
       info = -5;
     endif
   endif