## GNU Scientific Library - Bugs: bug #42219, Division by zero in...

You are not allowed to post comments on this tracker with your current authentification level.

## bug #42219: Division by zero in "gnewton" when "f" and "fdf" differ

 Submitted by: Patrick Alken Submitted on: Mon 28 Apr 2014 03:42:25 PM UTC Category: Runtime error Severity: 3 - Normal Operating System: Status: None Assigned to: None Open/Closed: Open Release:

Mon 28 Apr 2014 04:10:42 PM UTC, comment #1:

The problem seems to be in this line of gnewton_set() in multiroots/gnewton.c:

122 state->phi = enorm(f);

In the example program, the initial guess x is the solution vector, so the F_DF call provides a vector f = 0.

Then in the first gnewton_iterate() call, on line

177 double theta = phi1 / phi0;

we divide by 0.

Patrick Alken <psa>
Mon 28 Apr 2014 03:42:25 PM UTC, original submission:

Hello GSL team,

The "gnewton" multiroot solver has a possibility of dividing by zero if the
provided functions "f" and "fdf" do not produce bitwise-identical output
given the same "x". This results either in non-convergence with a result
of "-nan" or, when common floating point exceptions are trapped, a floating
point exception.

The division takes place on line 177 of "multiroots/gnewton.c". This line
will be executed with "phi0=0" if the first call to "fdf" lands precisely
on the root, but the subsequent call to "f" at the same "x" yields a
residual greater than zero. (The bug may be triggered in other
circumstances as well, but this is the case we encountered in our code.)

Why would "f" and "fdf" produce different results? The manual states that
"fdf" provides an optimization by computing both "f(x)" and "J(x)"
simultaneously. Some optimizations (such as subexpression elimination) can
lead to roundoff differences. More generally, with optimizing compilers
and new instruction sets, it is difficult to guarantee bitwise-identical
results. In this case, it was FMA instructions (both FMA4 and FMA3) that
triggered the bug, but we have also encountered reproducibility problems
related to AVX alignment. For the sake of example, here is the
optimization that yielded roundoff differences with FMA:

const double f0 = ww - c1_/(hh) - 1.0;

vs.

const double w2 = w*w;
const double h2 = h*h;
const double f0 = w2 - c1_/h2 - 1.0;

I have attached a small program (bug_gnewton.c) that triggers this bug,
based on the manual's example code for multiroot solvers. Currently, it
fails either with:

Type: gnewton
status = the iteration has not converged yet
root = -nan

or, if GSL_IEEE_MODE="trap-common", with:

GSL_IEEE_MODE="trap-common"
Type: gnewton
Floating point exception (core dumped)

The expected output is:

Type: gnewton
status = success
root = 1.000000e+00

I can work around this bug for now (by having "fdf" defer to "f" and "df"),
but it would be great if "gnewton" could be made more robust. We'd hate
for other clients to trigger this halfway into a multi-month simulation,
and unfortunately the other multiroot solvers have their own
division-by-zero problems (bug report coming soon). Thank you for looking
into this, and if you have any questions, just let me know.

- Curran

Patrick Alken <psa>

Attached Files
file #31265:  bug_gnewton.c added by psa (2kB - text/x-csrc)

Depends on the following items: None found

Items that depend on this one: None found

Carbon-Copy List
• -unavailable- added by psa (Submitted the item)
•

Please enter the title of George Orwell's famous dystopian book (it's a date):