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


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

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

Submitter:  Patrick Alken <psa>
Submitted:  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

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>
Group administrator
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 = w*w - c1_/(h*h) - 1.0;


    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:

    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>
Group administrator


(Note: upload size limit is set to 16384 kB, after insertion of the required escape characters.)

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


Depends on the following items: None found

Items that depend on this one: None found


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

    Follows 1 latest change.

    Date Changed by Updated Field Previous Value => Replaced by
    2014-04-28 psa Attached File- Added bug_gnewton.c, #31265

    Back to the top

    Powered by Savane 3.13-b921.
    Corresponding source code