[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
[Bug-gsl] [bug #42219] Division by zero in "gnewton" when "f" and "fdf"
From: |
Patrick Alken |
Subject: |
[Bug-gsl] [bug #42219] Division by zero in "gnewton" when "f" and "fdf" differ |
Date: |
Mon, 28 Apr 2014 15:42:26 +0000 |
User-agent: |
Mozilla/5.0 (X11; Linux i686 on x86_64; rv:28.0) Gecko/20100101 Firefox/28.0 |
URL:
<http://savannah.gnu.org/bugs/?42219>
Summary: Division by zero in "gnewton" when "f" and "fdf"
differ
Project: GNU Scientific Library
Submitted by: psa
Submitted on: Mon 28 Apr 2014 03:42:25 PM GMT
Category: Runtime error
Severity: 3 - Normal
Operating System:
Status: None
Assigned to: None
Open/Closed: Open
Release:
Discussion Lock: Any
_______________________________________________________
Details:
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 = w*w - c1_/(h*h) - 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
_______________________________________________________
File Attachments:
-------------------------------------------------------
Date: Mon 28 Apr 2014 03:42:25 PM GMT Name: bug_gnewton.c Size: 2kB By:
psa
<http://savannah.gnu.org/bugs/download.php?file_id=31265>
_______________________________________________________
Reply to this item at:
<http://savannah.gnu.org/bugs/?42219>
_______________________________________________
Message sent via/by Savannah
http://savannah.gnu.org/
- [Bug-gsl] [bug #42219] Division by zero in "gnewton" when "f" and "fdf" differ,
Patrick Alken <=